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Abstract 



We present a new Monte Carlo algorithm for simulating quantum spin systems which is able to 
suppress the negative sign problem. This algorithm has only a linear complexity in the lattice size 
used for the simulation. A general description and a rigorous proof of its correctness is given. Its 
efficiency is tested on a simple 2-dimensional fermionic model. For this model we show that our 
algorithm eliminates the sign problem. 
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1 Introduction 



The study of lattice quantum spin systems is a difficult problem. These systems are im- 
portant in condensed matter physics since they can describe strongly correlated electrons. 
At present, the domain of applicability of the existing analytical methods able to analyze 
them is still very restricted. To overcome this difficulty numerous groups have tried to 
study these models using computational techniques. 

Two major techniques are currently used to study quantum spin systems numerically. 
The first is the exact diagonalisation method and the second is quantum Monte Carlo (for a 
review see Q). Usually, both methods have a complexity which grows expo nentially with 
the size of the system under study, so that meaningful results are difficult to obtain. In 
particular, quantum Monte Carlo can suffer from the negative sign problem, which becomes 
exponentially serious on large lattices and at small tempera ture. 

In this work we develop a new algorithm which is able to suppress the sign problem in 
a quantum Monte Carlo simulation. A general description of this algorithm and a rigor- 
ous proof of its correctness are presented. This algorithm was briefly presented in 0]. It 
is based on a well-controlled approximation which has a linear complexity in the lattice size. 

We have applied it to a simple 2-dimensional models with fermions for testing its ef- 
ficiency. This model has a serious sign problem if simulated with conventional quantum 
Monte Carlo algorithms, but if simulated with our algorithm the sign problem is completely 
eliminated. The generality of this algorithm allow us to apply it to any model, opening 
new perspectives in the numerical study of quantum spin systems. 

This paper is organized as following: first we define the formalism and prove some 
lemmas (section 2 and 3) that we need for the description of the algorithm. The algorithm 
is then discussed and a proof of correctness is given (section 3). Finally in sect ion 4 the 
results of the example are presented. 

2 Formalism for quantum spin systems on a finite lattice 

We consider a d-dimensional finite lattice A C Z d . Associated to each site of the lattice we 
consider a particle with a finite number M of internal degrees of freedom. A Hilbert space 
7i x is associated with each site xsAof the lattice and is isomorphic to C . We choose an 
orthonormal basis {e%} a ^j with / = {1, M} in the Hilbert space TC X . The Hilbert space 
Ha over the lattice A is given by the tensor product 

H A = (g)H x (1) 

x<=A 
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The set J1a of all configurations oo\ in A is defined as the assignment of an element a x to 
each lattice site x £ A. The set {e WA } WAg n A of all vectors 

«-a = (g) 4. (2) 

forms an orthonormal basis of TCa- 



A Fock space is constructed to be able to incorporate the statistics of the particles 

fp(H A ) = P0^ (3) 

n>0 

where TL\ is the n-fold tensor product with itself, TL\ = C and P is the projection onto the 
symmetric or antisymmetric subspace for bosons or fermions, respectively. 

1" xeA 

PFermi(euj A ) = -y ^ ( 4 ) 

where 7r is a permutation of the order of the points (x,a) € A x I and sgn(ir) is +1 if 
the permutation ir is even —1 if it is odd. We choose an arbitrary order of the points in 
A denoted by x -< x' if x proceeds x' and we define the order of the points in A x 7 by 
(x, a) ^ {x' , cr') if a; -< x' or if x = x' , a < a' . 



The annihilation and creation operators on J-p(7i\) are defined as 

Cxa = PttxoP (5) 



and 



with 



C = Pa* xa P (6) 

a^(e£ ® ... ® e£) := v^(«£,e£)e£ ® ... e£ 

<4r(e£ ® ... ® e*») := Vn + le* ® e£ ® ... ® e£ (7) 

where (e^,e^) denotes the scalar product of the vectors and e^J. Furthermore, we have 
ow|0 >:= and a* a \0 >= e x a where |0 > denotes the zero particle state (vacuum). 



For bosons and fermions the operators defined through (||) and (^) satisfy the canonical 
commutation relations. For bosons we have 

[Cxa, Cx'a'] = [c xo ., C*/ CT /] = 

[c^,c* v ] = (<£,e£) (8) 
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and for fermions 



{c xa ,C x 'a'} = {c* CT ,<w} = 

{Cxa, <&*>} = ( 9 ) 



An orthonormal basis of J-p(7i\) is given by the vector 



1 CT 1 '" x k a k 



. 91=1 q k =l 

= (i4 ia X xiai -M^r^) ^ |o > (io) 

lll = l n XitTi- 

where the (...)-< denotes that subscripts (xiai) of the braced factors are ordered as discussed 
above. 



We define the algebra of observables A\ as the algebra of operators generated by the 
identity and all monomials of even degree in the creation and annihilation operators as- 
sociated with the lattice sites x £ A. Two algebras of observables A\ 1 and A\ 2 satisfy 
locality 

V4i G A Al ,A 2 G Aa 2 : [A Xl A 2 ] = if A x n A 2 = (11) 

and inclusion 

A Al C Aa 2 if Ai C A 2 (12) 

The trace of an operator A G A\ is defined by TrA = J2 V < > where {\v >} is an 

orthonormal basis of J-p{7i\) of the form ([lC|). 



The dynamics of the system is described by a finite range Hamiltonian 

H=J2$b (13) 

where $g is a selfadjoint operator belonging to Ab defined on a bond B C A. In our 
work we consider only periodic boundary conditions. The definition of a bond B is then 
adapted to our boundary conditions. If the basis vectors (|l0|) are eigenstates of H then 
the Hamiltonian describes a classical spin system, otherwise it describes a quantum spin 
system. The partition function of the statistical system is given by 

Z = Tre~ pH (14) 

where (3 is the inverse temperature. The expectation value of an observable A G Ax is 
obtained by 

< A >= ^-TrAe' 1311 (15) 

Zj 
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3 Monte Carlo of quantum spin systems 



3.1 The Trotter formula 



Classical spin systems can be simulated in d-dimensional lattice using Monte Carlo algo- 
rithms. The weight < v\e~@ H \v > contributing to the partition function can be directly 
evaluated since the basis vector (|l0|) is an eigenvector of the Hamiltonian. The evaluation 
of the weights is essential for any Monte Carlo algorithm. 



For a quantum spin system a Monte Carlo simulation is not straightforward because 
the basis vector (10) is not an eigenvector of the Hamiltonian any more and the evaluation 



of the exponential in the partition function becomes a very complex task. In fact, if one 
wants to perform any calculation directly using the d-dimensional lattice state \v >, the 
request in storage and computer time grows like the dimension of the Fock space, which is 
exponential in the volume |A| of the lattice. The on ly practical way to evaluate the weight 
< v\e~ p h \v > is to apply the Trotter formula in a (d+l)-dimensional lattice, where the 
extra dimension is the discretisation of the inverse temperature (we call it time) 

Z = lim Tr (e'^Y (16) 

where e = 1/T. For a simulation we have to restrict the interval X = [0,T] n Z to a finite 
set, approximating in this way the partition function. The systematic errors introduced by 
the approximation are of order e 2 . 

We consider a decomposition of the Hamiltonian (|l3|) into a sum of q operators 

H = z2Hi (17) 

i=l 

where each of them is a sum of commuting selfadjoint operators 

Hi = J2^k (18) 

From the locality of the algebra of operators we can construct this sum so that the operators 
<& l B satisfy 

[* l Bi ,<] = Oifs i n J B l ', = 

v 

This means that for each label i = 1, ...q we consider the sum J2b over disjoint bonds. 



We can define a classical (d+l)-dimensional spin system starting from the Trotter for- 
mula and inserting a projector 

1 (i,t) = E M*) >< (19) 

Vi(t) 
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for each time point in t £ X and each exponential e P eHi , Here the sum J2 Vi (t) g° es over an 
orthonormal basis of J~p(J~Lh)- The partition function is then approximated by 



Z ~ E < v 

»(0) 



£<«(o)i 



n n^i K°)>= 
.t=i \i=i /. 



|u(0) > 



« Lt=l \i=l 

Here the sum X^fo) S oes over an orthonormal basis of Tp(Tl\). The sum 



E : = E 

i>i(0),...,i> 9 (0),...,MT),...,VT) 



goes over an orthonormal basis of J-p(J~i\) where T = [0. q x T '"] Z and 



(20) 



(21) 



replicates the Fock space .Fp(Ha) at each slice r G T. Using the locality of the algebra of 
observables and the decomposition of the Hamiltonian (|Iq ) in commuting operators we can 
factorize the contributions in (|20|) in a product of terms localized over the bonds Bi 



Z ~ E w ( 7 



(22) 



where the weight of the vector \v >£ JFp (Wa) is given by the product of terms 



with, for i > 0, 



and 



and, for i = 0, 



n nn^(MM)) 

t=0 i=l Bj 



t' = t, i' = i + 1 if 1 < i < g 
i' = t + 1, i' = 1 otherwise 



w Bi (v,(i,0)) 



< V(U)\< 

1 



i |«i(l) > if i = 1 
otherwise 



(23) 

(24) 
(25) 

(26) 



The weight we^v, (i,t)) is a real number. We can decompose the weight wp^v, (i,t)) in 
its modulus and its sign. 



w Bi (v,(i,t)) = sgn(w Bi (v,(i,t))) x \w Bi (v, (i, t))\ 



(27) 



The evaluation of the weight \wBi(v, (i,t))\ requires to perform the exponential e s »' of 
the operator <& B f . Because this operator is local, it is a finite matrix of dimension equal to 
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the dimension of the Fock space localized on the bond Bi. The computational complexity of 
this operation is then very small. The sign of the weight sgn(wBi(v, (i,t))) is not generally 
a local quantity. For fermionic systems its evaluation requires the calculation of the sign 
of the permutation ir associated to the order in (|To|) of the < Vi(t)\ and |?v(t') > vectors. 
Also this operation is relatively easy. For bosonic systems the sign can also be negative but 
its evaluation remains local. 



3.2 Example 

As an example we consider a simple model of free quantum spin systems^. We consider 
fermions living on the sites of a spatially 2-dimensional square lattice with periodic bound- 
ary conditions. We consider a particle with only one internal degree of freedom I = {1}. 
We chose {e^}o- e / to be an orthonormal basis of the Hilbert space TC X . The Fock space is 
constructed as in (|3|). Creation and annihilation operators c* and c x anticommute 

{c*,c*} = 0, {c x ,c y } = 0, {c*,Cy} = S xy . (28) 

for x, y E A. We consider the Hamilton operator 

H = ^«c x + c* x+ ^ x+l - c* x c x+l - c* x+} c x ), (29) 

x,i 

where i is the unit vector in the i-direction. The model is trivial and can be solved in 
momentum space. However, when one tries to simulate it with a Monte Carlo algorithm it 
shows from the algorithmic point of view all the characteristic features of more complicated 
quantum spin systems. We can approximate the partition function Z with the Trotter 
formula ( |20|) by decomposing the Hamiltonian into four pieces H = H\ + Hi + H$ + H± 
(see eq. ([17|)) and 

H\= h x ,ii H 2 = 2J H%= 2J hx,i, H4 = h Xi 2, (30) 

x=(2m,n) x=(m,2n) i=(2m+l,n) x=(m,2n+l) 

where h Xj i = c*c x + c* x+ ~.c x+ - i — c* x c x+ ~ i — c* x+ -.c x (see eq. ([Tq)). The individual contributions 
to a given Hj commute with each other, but two different Hj do not commute. 



Using ( |20D we can write the grand canonical partition function 

Z = Tre-^ H -^ = UmT-,ooTr[e-^-i *) e -P<H 2 -1 iV) e -^(^ 3 -f N) e -fk(H<->t N) ]T 

(31) 

where [i is the chemical potential and = J2xgA n x with n x = c x c x is the particle number 
operator. Inserting a complete set of Fock states between the factors and using the locality 
of h x ,i for bonds < x,x + i > we obtain matrix elements of the exponential (p4|). 

e -/?e(Vi-£(n„+n c+ ,)) = e -^e(l-f) x 

/ exp(^e(l - f )) \ 

cosh(/3e) sinh(/3e)^ 

sinh(/3e)SI cosh(/3e) 

V exp(-/3e(l - f )) J 



(32) 



2 This example was analyzed by quantum Monte Carlo in 
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where the 4x4 matrix is in the Fock space basis 1 00 >, |01 >, 1 10 >, |11 > defined on 
the bond < x, x + i >. 51 is the sign of the weight which has to be evaluated looking at 
the permutations of the fermions with respect to the given order of the basis vectors (|l~0|). 
The partition function is approximated by a classical spin system over occupation numb ers 
n(x, t) =0,1 (here r labels the T ■ q time slices in T with q = 4) with periodic boundary 
conditions. The systematic errors due to the finite T are of order e 2 . 

^ — II a2 \w(n)\ sgn(w(n)). (33) 
z,t n(a,T)=0,l 

The modulus factor takes the form 

l u '( n )l = II w[n(x, r), n(x + 1, r), n(x, r + 1), n(x + 1, r + 1)] 

x=(2m,n),T=Ap 

x JJ w[n(x, r), n(x + 2, r), n(x, r + 1), n(x + 2, r + 1)] 

x=(m,2n),T=4p+l 

x Yl w[n(x, t), n{x + 1, r), n(x, r + 1), n(x + 1, r + 1)] 

z=(2m+l,n),T=4p+2 

x ]~J w[n(x,r),n(x + 2,r),n(x,r + 1), n(i + 2, r + 1)], 

x=(m,2n+l),r=4p+3 

(34) 

with w[0, 0, 0, 0] = 1, w[l, 1, 1, 1] = e~ 2/3e(1 - 4 ), w [0, 1, 0, 1] = w[l, 0, 1, 0] = e _/3e(1 - f } cosh(/3e), 
w[0, 1,1,0] = io [1, 0,0,1] = e~ /3t ( 1 ~ 4) sinh(/3e). All other values are zero. The occupation 
numbers n(x, r) = 0, 1 are variables interacting with each other via the time-like plaquette 
couplings w[n(x, r), n(x + i, r), n(x, r + 1), n(x + i, r + 1)]. 

Each state is weighted by a sign factor which arises from the fermionic statistics. 
The sign factor sgn[w(n)} is a product of terms sgn[n(x, r),n(x + i, r), n(x, r + 1), n(x + 
i, r + 1)] associated with each plaquette interaction. One has sgn[Q, 0, 0, 0] = sgn[l, 1, 1, 1] 
= sgn[0, 1,0, 1] = sgn[l, 0, 1, 0] = 1. A nontrivial sign ±1 may arise only for plaquette 
interactions of type [0, 1, 1, 0] and [1, 0, 0, 1]. 



3.3 The sign problem 

The decomposition of the Hamiltonian in local terms ( |i~7| , |l8| ) allows us to treat the new 
(d+l)-dimensional spin system as a classical spin system with state space J-p(7i\), par- 
tition function (^) and weight (22). The evaluation of the weight is of small complexity, 
contrary to the original d-dimensional quantum spin system (|l~4T). The price to pay is that 
the new classical system has a partition function (22) which has not generally positive 
semidefinite weights w(v). This fundamental difficulty is usually referred to as the "neg- 
ative sign" problem. It is not related to any approximations in the Monte Carlo scheme 
but it describes the fact that the statistical error of the observables can become very large, 
increasing exponentially in the inverse temperature (3 and lattice volume |A x T\. 
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Definition 1 A classical observable is an operator A £ A-a acting on the Fock space 
Fp(1~La) defined at t = and diagonal with respect to the basis of the form (|To|). We 
denote by A(v) the matrix element A(v) =< v(0)\A\v(0) >. 

For simplicity we restrict our discussion to classical observables. With some minor 
modifications our algorithm is applicable also for observables which are not diagonal with 
respect to the basis (|To|). We emphasize that, usually, the interesting obs ervables are 
calssical. 



The expectation value of a classical observable can be written 

<A>= Y. v A{v)w{v) = EyAjv) \w(v)\sgn(w(v)) 
Y,v w i v ) T,v \w(v)\sgn(w(v)) 

A Monte Carlo algorithm needs a positive semidefinite weight to be able to construct a 
Markov process. Redefining the observable incorporating the sign in it A(v) = sgn(w(v))A(v) 
we can simulate a positive semidefinite classical spin system with statistical weight |u)(u)| 
and correct the measurement of the observable 

< A >„= S^WkWI x 1 =; <^> M (36) 

z2 v FWl < sgn{v) >\ w \ < sgn >\ w \ 

where < ... > w and < ... >\ w \ denote the averages taken with the weights w and 
respectively. If the average sign is small there will be large cancelations making an accurate 
evaluation of < A > impossible. 



3.4 Equivalent Monte Carlo algorithms 

We consider a Monte Carlo algorithm cj)\ w \ which produces a Markov chain with state space 
{\v >G J-p(H.A)} and equilibrium distribution We assume that this algorithm 

generates a state \v' > from a state \v > with a transition probability T(v' <— v) 

cj) H : \v >^ \v' > (37) 

We assume stationarity 

Y,T(v' ^ v)\w{v)\ = \w(v')\ (38) 

V 

and irreducibility (ergodicityQ) 

T(v' ^v)>0, y\v >, \v' >G Fp(Hh) with |to(t;)| > and \w{v')\ > (39) 

to ensure the convergence to equilibrium. This algorithm can be realized in different ways. 
The explicit realization of it is for the moment irrelevant to the discussion. We denote the 

3 Notice that only on Markov chains with finite state space, "ergodic" is used in the physics literature as 
a synonym for "irreducible" (see |H|, Section 2.4). On Markov chains with general state spa ce they have a 
different meaning (B, p. 169). 
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expectation value measured with the algorithm cj)\ w i by eq. (|36|). 

If the average of the sign is small this algorithm will suffer from the sign problem as 
discussed above. It is important to notice that the value A(v) of a classical observable A 
does not depend on the complete state \v > but only on the part living on the time slice 
at t = 0. This is a crucial property that we will use for constructing an algorithm whi ch 
is able to reduce the sign problem by substantially increasing the average sign. We first 
introduce some simple definitions and lemmas which allow us to construct this algorithm. 

Definition 2 A mapping 

g : Tf(H A ) - Fp(H a ) 

\v >h-> g\v > (40) 

is called observable preserving if for any classical observable A its value A(v) is g invariant 
A(v) = A(gv). Such a mapping is easily constructed: we can require that it does not change 
the state \v > at t = 0. 



Lemma 1 We consider a set Q of observable preserving and bijective mappings g 
define the new weight 

w(v,g) = ~(w(v) +w(gv)) 
then we obtain for the expectation value (^) of a classical observable 

EyA{v)w{v) = Eggg EvMv)w(v,g) 

Eutt'(u) E 9e gE„*KsO 

Proof : Using the bijective property of g, it is easy to see that after summing over all states \v ><G 
J-p(Ha) we obtain J2 V w ( v ) = Et, W (9 V )- We denote by \Q\ the number of elements in Q and after 
summation over all mappings g £ Q we obtain that ~^2 V w(v) = A E 9 eg E« w(v, g). Using the fact 
that g is observable preserving (that means A(v) = A(gv)), in analogy as before we can see that 
Et, A(v)w(v) = ^ E 9ee E, A(v)w(v,g) so that © holds. □ 

Definition 3 We call two Monte Carlo algorithms equivalent if for all classical observables 
A their expectation values measured with both algorithms are equal. We call two Monte 
Carlo algorithms 5-quasi equivalent if for all classical observables A their expectation values 
measured with both algorithms are equal up to systematic errors of order 5. 

Using this last lemma || it is clear that the expectation value of A measured with a Monte 
Carlo algorithm (j>i w i is the same as the one measured with a Monte Carlo algorithm so 
that and are equivalent. The set of mappings Q can contain an arbitrary number 
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We 
(41) 

(42) 



of elements. We suppose that we can construct a mapping g so that the average of the sign 
if measured with (f>^ and w defined as in (|4~i"|) satisfies 

< sgn(v) >|^| > < sgn(v) >\ w \ (43) 

where < ... >\w\ and < ... >\ w \ mean the expectation taken from the \w\ and distribu- 
tions, respectively. In the case that the average sign substantially increases then it is evident 
that the sign problem is suppressed. Notice that the sign is not a classical observable be- 
cause it is a global quantity depending on the state over the complete lattice A x T . Two 
equivalent Monte Carlo algorithms can thus have different sign expectation values. The 
sign prob lem is now shifted to the search of a mapping with the desired property. Such a 



mapping can be constructed looking for clusters of spins to be flipped so that (43) is satisfied 
with high probability during the Monte Carlo process. This is a non trivial problem, since to 
ensure the bijectivity of the mapping the required work grows exponentially in the volume 
of the lattice AxT because a deterministic cluster search is needed. Apart from some triv- 
ial example, where the amount of work is not too excessive, this approach is not convenient. 



However, the problem can be solved by constructing a mapping which is bijective with 
probability 1 — 5 at each Monte Carlo iteration where 5 is very small. We will see that the 
construction of such a mapping has a complexity which is linear in the volume of the lattice 
because a randomized cluster search can be used in connection with a hashing technique for 
testing the bijectivity. Since the mapping is not always bijective, during the Monte Carlo 
process a systematic error is produced, so tha t is only 5-quasi equivalent to (j)\ w \. If 
we choose 5 to be smaller than the statistical error and the systematic error e 2 introduced 
by the Trotter formula, the precision of our algorithm is sufficient. In the next subsection 
we will show how this mapping can be constructed explicitly. 

3.5 Reduction of the sign problem 

We want to construct a mapping g which has the property (^3|) and can be used to construct 
a Monte Carlo algorithm J-quasi equivalent to 4>\w\- First we introduce some concepts which 
allow us to construct this mapping. 

Definition 4 A mapping 

g : ^(H A ) - ^(H A ) 

\v >^ g\v > (44) 

is called compatible in \v > if w(y) ^ then w(gv) ^ 0. 

Definition 5 A state \v > of the Fock space Jp(Ji.x) can be rewritten as 

\v>=: (g) \n a (x,r)> (45) 

(x,t)gAxT 
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where a £ I is defined in section 2. A flip e,( x ,t) w hh { x i T ) £ A x T is a mapping from 
Fp{T~Ca) to J r p(Ti.A) which locally transforms |n CT (x,r) > in some \n' a ,(x,r) > so that 
S/ X)T ) o S( XjT ) = 1. The composition of flips over a subset Q C A x T is denoted by Sq. A 
subset C C A x T is called a cluster of |t> > if He is compatible in |u >. A subset C C A x T 
is called a preserving cluster of \v > if He is compatible in \v > a nd observable preserving. 
We call macros of \v > the set of preserving clusters of a state \v > with respect to a flip 
E and denote it by M(v,E). For a given flip we just denote it by M(v) if the context 
allows that. We defin e an arbitrary order of clusters in A4(v) and denote it by C -< C if C 
precedes C. We consider the empty cluster 6 M(v) to be ^ C, VC G A4(f). 

We now proceed to the construction of the mapping g with the desired properties. 

Algorithm 1 We consider a hashing function h which assigns in an arbitrary way a non 
vanishing integer label to any state 

h : J%(Ha) - 2 C Z 

\v >^ h(v) (46) 

We suppose that we can store a hashing table H table with \Z\ integer entries. Here \Z\ 
denotes the number of different labels and \Z\ < dim(J-^> (7i\)). Given a state \v > we 
define g\v > with the following procedure: Let M! C M.{v) be an arbitrary subset of M.{v) 
containing the empty cluster 0, then 

select the first C G M' 

repeat 

if (condition=true) then 

g\v >= Ec\v > 
otherwise 

g\v >= \v > 

select next C £ M' 
end if 

until (condition=true) or (C = 0) 

where we define the 

condition = {(w(v) +w(E c v) > 0) and 0(v,E c v)} (47) 

The boolean function 0(v,v f ) is defined by the following procedure 

0(v,v') : if H table {h(v')) = then H taUe {h{v')) := h(v) endif 
if H ta bie{h(v')) / h(v) then output 0=false 
else output C=true endif 
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The selection of the clusters in the macros follows the chosen order of the clusters. Practi- 
cally there is no need to find all the clusters of a macros. One can select a point (x, r) G AxT 
and construct a cluster starting from it. Duri ng the construction a fixed list of random 
numbers can be used. It is important, however, that this list remains always the same every 
time one applies this procedure to a state \v >. Changing the list of points or the list of 
random numbers is equivalen t to select a new mapping g in Q. If the constructed cluster 
does not satisfy the condition (47) the next point in A x T can be selected and a new cluster 
constructed and this search is repeated until the conditio n (47) is satisfied. If the condition 
(|47| ) is never satisfied the procedure can be stopped, for example, when all the points in 
AxT are tested once, and the original state is returned as the result. In this way the s 
earch of the cluster has a linear average complexity in the size of the lattice. 



Theorem 1 Let Q a set of mappings defined by the algorithm |T[ A weight w is defined by 

w{v, g) = ~(w(v) + w(gv)) (48) 

We consider a Monte Carlo algorithm cft^ which produces a Markov chain with state space 
{ip = (\v >,g) € T"p(Wa) x G} and equilibrium distribution \w(v, g)\. We assume that this 
algorithm generates a state ip' = (\v' >,g') from a state ip = (\v >,g) with a transition 
probability T{ip' <— ip). We assume stationarity 

Y, T W «- ^)\w{v,g)\ = \w{v',g')\ (49) 

and irreducibility 

T(ip' «- ip) > 0, \/ip,ip' 6^p(Ha) X G with \w(y,g)\ > and \w(y',g')\ > (50) 
to ensure the convergence to equilibrium. 

Then this Monte Carlo is 5-quasi equivalent to <fi\ w \ with 5 = and 

< sgn(v) >\yj\ > < sgn(v) >\ w \ (51) 



Proof : The flip is defined observable preserving (see definition |5|) so that all g are also. The 
boolean function O in the condition ( |47| ) guarantees us that the mapping g maps a state \v > 

bijectively to a state \v' > with probability 1 — (j^J ■ Because of lemma |] it is then clear that 

4>\w\ is equivalent to cf)\ w \ up to errors of order in the average of an observable, so that the two 

algorithms are (5-quasi equivalent with <5 = "T^T " ^ e conc '' 1: ' on (0) m algorithm ^[guarantees us that 

{sgn(w(v) + w(Ec(v))) = 1 if a cluster C satisfying "1 
condition is found in M' , > > sgn{w(v)) 
sgn(w(v)) otherwise J 

so that after average (^J) is satisfied. □ 

It is important to notice that the systematic error of order 5 introduced by this algorithm 
can be tuned by increasing the dimension \Z\ of the hashing table. 
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3.6 Monte Carlo algorithm for the weight w. 

A Monte Carlo algorithm for the spin system described by the weight w(v,g) needs a 
method for updating the state tp to a new state tp' with a transition probability T(tp' <— tp) 
which satisfies stationarity (??) and irreducibility (|50|), The modulus of the weight \w\ is 
not local, contrary to (see section 3.1), because after adding the term w(3cv) in w 

we can not factorize w anymore in a product like (^3|). To find an efficient updating method 
applicable directly to w is a difficult task. One can, however, update the states \v > using 
the Monte Carlo algorithm <p\ w \ which is supposed to be known and then correct the weight 
distribution using an accept/reject global Metropolis test. 

Algorithm 2 Suppose we know an Monte Carlo algorithm <j)\ w \ for the distribution 

For simplicity we suppose that the Monte Carlo algorithm <f>i w i satisfies irreducibility and 

the detailed balance condition 

T(v' <- v)\w(v)\ = T{v <- v')\w(v')\ (52) 

where T(y' <— v) is the transition probability of <p>\ w \. It is clear that stationarity follows 
from this condition. We consider a set of mappings Q constructed as in algorithm |l[ We 
define the Monte Carlo algorithm <f>\$\ for the equilibrium distribution \w(v,g)\ and the 
state space {tp = (\v >,g) € Jp^Hh) x G} by the procedure 



input the state tp = (\v >,g) 
generate \v' >= <fi\ w \\v > using (p\ w t 
select randomly a mapping g' £ Q 

\w(v',g')\ ■ \w(v)\ 



PA(tp',tp) = min fl, 



\&{v,g)\ ■ \w(v')\ 

accept tp' = (\v' >, g') as the new state with probability P^itp' , tp) 



Theorem 2 Let ^i^i be defined as in algorithm ^. Then (pi^l satisfies stationarity ( |49| ) and 
irreducibility (|50|). 



Proof : Stationarity can be proven by showing that the detailed balance condition 

f (V/ «- 1>)\wU>)\ = f{il> i- 1>')\w{i/>')\ 

is satisfied by <f>^ where 

TU>'^1,)=T(v'<-v)-Pa(iI>',iI>) 

is the transition probability for 4>\w\- We suppose that Pa{"4>' \tp) < 1 for the states \tp > and \tp' >. 
Using the property that Pa{^, tp') = 1 if Pa{"P\ tp) < 1 and eq. ( p^ ) we have 

f(tp' - tP)\w(tP)\ = T(v> <- v) ■ P A (tP',tP)\w(tP)\ = T(v> - ^ fcj^i'/ffjl NWl = 
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= T(v <- w ) u,(t;) 1— — = T(u <— w ) \w(v ) , , ,.. = 
\w{v')\ \w(v ')\ 

= T(v «- */) W)l = r (« <- «') ' /MV^O W)l - 

= t(v^v')HV>')I 

The case Pa(iP',^) = 1 and Pa{iP,iP') < 1 is analogous. Irreducibility is clear because w(ip) ^ 
when w(v) ^ so that Pa (?/>', "0) > and because <f)i w i satisfies irreducibility. □ 



If the set of mappings Q contains only one element then the dimension of the hashing table 
has to be larger than the number of Monte Carlo iterations one desires to perform. In this 
way one avoids too many collisions in the hashing table. This, of c ourse, is doable, but 
requires a lot of memory. If, however, the set of mappings Q contains a huge amount of 
elements then the hashing table can be small because the algorithm uses a selected mapping 
only for a short time and then selects a new one according to the acceptance probability 
Pa- If the set Q is big enough, the probability that the algorithm selects the same mapping 
twice is infinitesimal so that there is no need to store the history of the hashing tables of 
old mappings. 



Theorems |l] and [2] show that a Monte Carlo simulation can be performed using and 
a classical observable A can be measured using 

<A>m= i^2M + o(-L) (53) 
< sgn >| t5 | \\Z\J 

If the condition ( p7|) is satisfied with high probability during the Monte Carlo process then 
the negative sign problem is eliminated. The systematic error O ( produced by the 
hashing technique used for checking t he bijectivity of the mappings g can be explicitly 
measured durung the simulation. 



4 Application to the example 

We apply our algorithm |2| to the example presented in section 3.2. We use for the updating 
of the states \v > with respect to the distribution |itf(i>)| a standard loop algorithm M, |j 
which we denote as 4>\ w \. A description of this algorithm applied to the example of sectio 
n 3.2 is given in 1| . The reduction of the sign problem is then realized using our algorithm 0. 

The loop algorithm is essentialy a cluster algorithm where a cluster C in the macros 
Ai(v) is selected with a certain probability so that an update of \v > realized by a flip 
Sc satisfies detailed balance. For c ompleteness, we briefly describe the loop algorithm we 
have used in the example. We define the flip ^(x,t) m our example so that the occupation 
numbers n(x,r) of points on the cluster are changed from to 1 and vice versa: 

3(x,t)|1 >= 1° > 
h (*,t)|0 >= |1 > 
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The clusters in M{v) are constructed searching closed loops with the following algorithm. 



Algorithm 3 This algorithm finds a cluster C in A4(v): 

1. To start a loop one first selects a lattice point (x, r). 

2. The occupation number n(x, r) participates in two plaquettes, one before and one 
after r. For n(x, r) = 1 we consider the plaquette at the later time and for n(x, r) = 
we consider the plaquette at the earlier time. 

3. The corresponding plaquette is characterized by the occupation numbers of four points 
in A x T. One of these points will be the next point on the loop. 

• For a plaquette [0,0,0,0] or [1,1,1,1] the next point is with probability p\ the 
time- like nearest neighbor of (x,t), and with probability 1 — p\ the next-to- 
nearest (diagonal) neighbor of (x, r) on the plaquette. 

• For a plaquette [0, 1,0, 1] or [1,0, 1,0] the next point on the loop is with proba- 
bility p2 the time-like nearest neighbor, and with probability l—p2 the space-like 
nearest neighbor of (x,t). 

• For a plaquette [0, 1, 1,0] or [1,0,0, 1] the next point is with probability p% the 
diagonal neighbor, and with probability l—ps the space-like nearest neighbor of 
(x,r). 

4. Once the next point on the loop is determined the process is repeated from 2 until 
the loop closes. 

5. The points on the closed loop determine the cluster C. 

The sets of all clusters determined by this algorithm for some arbitrary values of Pi,P2,P3 
is only a subset of the macros M{v). This subset is, however, sufficient for the construction 
of our Monte Carlo algorithm. 

If the start point in algorithm ^ is chosen randomly and the probabilities pi,P2 and P3 
are chosen in the following way 

(Pl\ ( 1(1 + e"*) \ 

p 2 \ = \ i(l + e-^)/cosh(/3e) (54) 
\P3 ) \ (l-i(l + e-^))/sinh(/3e) ) 

an update of \v > realized by a flip He obeys detailed balance ||. The part of the weights 
proportional to exp(— /3e(l — ^)) is taken into account by a global Metropolis step. In the 
Metropolis step the cluster is flipped with probability p = min(l, exp(/3(4 — n)W(C))) where 
W(C) = ■^J2(x,T)£c(^ n ( x ^ T ) ~ -0- This defines the loop algorithm <fi\ w \. 
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The search method of the clusters in algorithm |l| can be performed using algorithm || 
choosing the starting point (cc, r) from a list of points in A x T and a list of random num- 
bers. There, the choice of the probabilities pi,P2 and p% is arbitrary. However, to obtain 
a good acceptance < Pa > it is convenient to use the probabilities defined in (|54|). The 
weight w(v) is defined by algorithm [l] and eq. (|48|). For our runs we have used a hashing 
table with 10000 entries. The loop algo rithm <j)\ w \ applied in algorithm || defines our Monte 
Carlo algorithm (j)^ . In our simulation we have not found any collision in the hashing table 
so that the systematic error 5 on the observables is zero. 

The model in our example is trivial and can be solved in momentum space. This allows 
us to test our algorithm. By introducing c* = J2 X eyi v{w x ) c xi c p = x^x ex ^P(~^P x ) c xi 
the Hamiltonian in momentum space becomes H = J2pP 2c p c p with pi = 2sin(j*j/2). In the 
grand canonical ensemble the expectation value of the occupation number is given by 

{«.) = >K eM-W - „jv))] = ^ £ TT1 ^1 W - W) (55) 

In Table 1 we present the results of the Monte Carlo simulations performed with the 
loop algorithm and with our algorithm § and we compare the obtained results with the 
exact solution (^). We have applied the algorithms for various values of (5 and [i at fixed 
lattice spacing (3e = 1/16. The results of both algorithms agree with the exact results within 
the error bars. It is evident that the sign problem becomes severe for the loop algorithm 
when the temperature is lowered or the chemical potential is increased. However, the sign 
remains always positive for our algorithm. 



5 Conclusion 

We have presented a new algorithm for significantly improving quantum Monte Carlo sim- 
ulations of models plagued by the negative sign problem. Its complexity is only linear in 
the volume of the lattice used for the simulation. The generality of this algorithm allows 
us to apply it to any quantum spin system. The efficiency of this algorithm was tested on 
a simple fermionic model. In this example the sign problem is solved. A more exhaustive 
analysis of the dynamics of the algorithm is under study. Applications of it to more physi- 
cally interesting models are also planned. 
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p 




Vi A i 


4T 


<n>o 


< n > L 


< n > E 


< sgn >o 


< sgn > L 


< Pa > 


0.5 


2 


8 


32 


0.304(2) 


0.305(3) 


0.3049 


1.0(0) 


0.690(4) 


0.64(1) 


0.5 


4 


8 


32 


0.501(2) 


0.498(3) 


0.5000 


1.0(0) 


0.591(4) 


0.64(1) 


1.0 


2 


8 


64 


0.233(2) 


0.230(40) 


0.2321 


1.0(0) 


0.048(6) 


0.612(9) 


1.0 


3 


8 


64 


0.357(2) 


? 


0.3568 


1.0(0) 


-0.003(6) 


0.60(1) 


1.0 


4 


8 


64 


0.499(2) 


? 


0.5000 


1.0(0) 


0.002(8) 


0.59(1) 


2.0 


2 


8 


128 


0.193(2) 


? 


0.1956 


1.0(0) 


0.004(9) 


0.58(1) 


2.0 


4 


8 


128 


0.498(2) 


? 


0.5000 


1.0(0) 


-0.001(9) 


0.58(1) 


1.0 





12 


64 


0.067(2) 


0.069(2) 


0.0685 


1.0(0) 


0.703(4) 


0.64(1) 


1.0 


1 


12 


64 


0.134(2) 


0.139(5) 


0.1351 


1.0(0) 


0.295(8) 


0.61(1) 


1.0 


2 


12 


64 


0.232(2) 


0.20(10) 


0.2321 


1.0(0) 


0.018(6) 


0.60(1) 


1.0 


3 


12 


64 


0.358(2) 


? 


0.3568 


1.0(0) 


0.001(3) 


0.59(1) 


1.0 


4 


12 


64 


0.501(2) 


7 


0.5000 


1.0(0) 


-0.001(4) 


0.58(1) 


1.0 


5 


12 


64 


0.642(2) 


7 


0.6432 


1.0(0) 


-0.002(5) 


0.59(1) 


1.0 


6 


12 


64 


0.768(2) 


7 


0.7679 


1.0(0) 


0.001(4) 


0.58(1) 


1.0 


7 


12 


64 


0.863(2) 


7 


0.8649 


1.0(0) 


-0.001(4) 


0.55(1) 



Table 1: Results from the Monte Carlo simulations for different parameters. The simulations are 
done with the loop algorithm 4>\w\ an d our algorithm [|. The expectation values of the number 
operator n and the sign obtained with the loop algorithm are denoted by < n >l and < sgn > i , 
the ones obtained with our algorithm |^ are denoted by < n >o and < sgn >o- The mark "?" 
indicates that the measurement of the expectation value was impossible, because the statistical 
error exceeds the expectation value. The exact values of < n > are denoted by < n >e- The 
average of the acceptance < Pa > of the algorithm || is also measured. 
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